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Abstract 

We  review  the  development  of  diffuse-interface  models  of  hydrodynamics  and  their 
application  to  a wide  variety  of  interfacial  phenomena.  These  models  have  been  ap- 
plied successfully  to  situations  in  which  the  physical  phenomena  of  interest  have  as- 
sociated with  them  a length  scale  commensurate  with  the  thickness  of  the  interfacial 
region,  (e.g.  near-critical  interfacial  phenomena  or  small  scale  flows  such  as  those  occur- 
ring near  contact  lines),  and  fluid  flows  involving  large  interface  deformations  and/or 
topological  changes  (e.g.  breakup  and  coalescence  events  associated  with  fluid  jets, 
droplets,  and  large- deformation  waves).  We  discuss  the  issues  involved  in  formulating 
diffuse-interface  models  for  single- component  and  binary  fluids.  Recent  applications 
and  computations  using  these  models  are  discussed  in  each  case.  Further,  we  address 
issues  including  sharp-interface  analyses  that  relate  these  models  to  the  classical  free- 
boundary problem,  computational  approaches  to  describe  interfacial  phenomena,  and 
models  of  fully-miscible  fluids. 
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1.  INTRODUCTION 

The  nature  of  the  interface  between  two  fluids  has  been  the  subject  of  extensive  investigation 
for  over  two  centuries.  Young,  Laplace  and  Gauss,  in  the  early  part  of  the  1800’s,  considered 
the  interface  between  two  fluids  to  be  represented  as  a surface  of  zero  thickness  endowed 
with  physical  properties  such  as  surface  tension.  In  these  investigations,  which  were  based 
on  static  or  mechanical  equilibrium  arguments,  it  was  assumed  that  physical  quantities, 
such  as  density,  were,  in  general,  discontinuous  across  the  interface.  Physical  processes  such 
as  capillarity  occurring  at  the  interface  were  represented  by  boundary  conditions  imposed 
there  (e.g.  Young’s  equation  for  the  equilibrium  contact  angle  or  the  Young-Laplace  equation 
relating  the  jump  in  pressure  across  an  interface  to  surface  tension  times  curvature).  Poisson 
(1831),  Maxwell  (1876)  and  Gibbs  (1876)  recognized  that  the  interface  actually  represented 
a rapid  but  smooth  transition  of  physical  quantities  between  the  bulk  fluid  values.  Gibbs 
introduced  the  notion  of  a dividing  surface  (a  ‘surface  of  discontinuity’)  and  surface  excess 
quantities  in  order  to  develop  the  equilibrium  thermodynamics  of  interfaces.  The  idea  that 
the  interface  has  a non-zero  thickness  (i.e.  it  is  diffuse)  was  further  developed  in  detail  by 
Rayleigh  (1892),  and  van  der  Waals  (1893)  who  proposed  gradient  theories  for  the  interface 
based  on  thermodynamic  principles.  In  particular,  van  der  Waals  gave  a theory  of  the 
interface  based  on  his  equation  of  state  and  used  it  to  predict  the  thickness  of  the  interface 
which  he  showed  became  infinite  as  the  critical  temperature  is  approached.  Korteweg  (1901) 
built  on  these  ideas  to  propose  a constitutive  law  for  the  capillary  stress  tensor  in  terms  of  the 
density  and  its  spatial  gradients.  These  original  ideas  have  been  developed  further  and  refined 
over  the  past  century  and  we  turn  the  reader  to  the  work  of  Rowlinson  (1979)  and  Rowlinson 
and  Widom  (1989),  who  provide  thorough  discussions  of  the  historical  perspectives  and 
complete  references  to  the  early  work  on  interfacial  and  capillary  phenomena. 

The  notion  of  a diffuse  interface  and  the  use  of  a capillary  stress  tensor  to  model  the 
interface  between  two  fluids  and  the  forces  associated  with  it  are  of  central  importance  to 
the  topics  under  consideration  in  the  present  paper.  Our  focus  here  is  on  the  use  of  diffuse- 
interface  models  which  fully  couple  these  notions  into  a hydrodynamic  description.  Such 
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models  have  been  used  to  understand  physical  and  hydrodynamic  phenomena  that  occur 
near  a fluid’s  critical  point.  Additionally,  developments  in  modern  computing  technology 
have  stimulated  a recent  resurgence  in  the  use  of  the  diffuse-interface  models  for  the  com- 
putation of  flows  associated  with  complex  interface  morphologies  and  topological  changes, 
such  as  droplet  breakup  and  coalescence  and  the  highly  nonlinear  development  of  classical 
hydrodynamic  instabilities. 

In  the  classical  fluid  mechanical  approach,  the  interface  between  two  immiscible  fluids  is 
modeled  as  a free  boundary  which  evolves  in  time.  The  equations  of  motion  which  hold  in 
each  fluid  are  supplemented  by  boundary  conditions  at  the  free  surface  which  involve  the 
physical  properties  of  the  interface.  This  formulation  results  in  a free-boundary  problem 
(Lamb  1932,  Batchelor  1967,  Lighthill  1978,  Drazin  & Reid  1981,  Davis  1983). 

Specifically,  in  the  free-boundary  formulation  it  is  assumed  that  the  interface  has  associ- 
ated with  it  a surface  tension,  which  on  applying  a stress  balance  at  the  interface  gives  rise 
to  the  interfacial  boundary  condition 

a-  • n|+  = 7/Cn,  (1) 

which  relates  the  jump  in  the  stress  across  the  interface  to  the  interfacial  curvature.  Here 
cr  is  the  stress  tensor,  n is  the  unit  vector  normal  to  the  interface,  7 is  the  surface  tension 
(here  assumed  to  be  constant)  and  JC  is  the  appropriately  signed  curvature.  In  addition,  an 
interface  between  two  immiscible  fluids  is  impermeable  in  which  case  conservation  of  mass 
across  the  interface  gives  that 

v ■h\_  = v-n\+  = Vn,  (2) 

where  v represents  the  velocity  of  the  fluid  and  Vn  is  the  normal  velocity  of  the  interface. 
Finally,  for  viscous  fluids,  there  is  continuity  of  tangential  velocity  across  the  interface 

[v  — (y  ■ n)n\\+  = 0.  (3) 

The  free-boundary  description  has  been  a very  successful  model  in  a wide  range  of  sit- 
uations, however  there  are  also  important  instances  where  it  breaks  down.  In  short,  as  a 
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physical  model  it  breaks  down  when  the  interfacial  thickness  is  comparable  to  the  length  scale 
of  the  phenomena  being  examined.  For  example,  (i)  in  a near-critical  fluid  the  thickness  of 
the  interface  diverges  at  the  critical  point  (Stanley  1971)  and  consequently  the  representation 
of  the  interface  as  a boundary  of  zero  thickness  may  no  longer  be  appropriate,  (ii)  the  motion 
of  a contact  line  along  a solid  surface  involves  a detailed  consideration  of  the  fluid  motion  in 
the  vicinity  of  the  contact  line,  and  may  require  the  treatment  of  length  scales  comparable  to 
that  of  the  interface  thickness,  and  (iii)  the  free-boundary  description  may  not  be  adequate 
for  situations  involving  changes  in  the  topology  of  the  interface  (e.g.  the  breakup  of  a liquid 
droplet),  since  these  processes  fundamentally  involve  physical  mechanisms  acting  on  length 
scales  comparable  to  the  interface  thickness.  In  addition  to  the  above  situations,  another 
difficulty  associated  with  the  free-boundary  formulation  arises  in  its  use  in  computational 
settings  when  the  free  boundary  shape  becomes  complicated  or  self-intersecting. 

Diffuse-interface  models  provide  an  alternative  description  in  the  face  of  these  difficulties. 
Quantities  which  in  the  free-boundary  formulation  are  localized  to  the  interfacial  surface  are 
recognized  to  be  distributed  throughout  the  interfacial  region.  For  example,  surface  tension 
in  the  classical  model  is  a representation  of  a distributed  stress  within  the  interfacial  region. 
In  this  spirit  a continuum  theory  of  the  interface  may  be  developed  where  the  reversible  part 
of  the  stress  tensor,  C,  that  is  associated  with  surface  tension  is  expressed  in  its  simplest 
form  as 

Coc  (/>V2/>  + l|W>|2)  I-VP®VP,  (4) 

where  p is  the  fluid  density  and  the  components  of  the  outer  product  Vp  0 Vp  are  given  by 
(dp/ dxi)(dp/ dxj).  The  existence  of  such  a stress  tensor,  often  called  the  capillary  tensor, 
was  first  described  by  Korteweg  (1901).  The  derivatives  of  the  density  that  appear  in  the 
stress  tensor  arise  from  the  non-local  interaction  of  the  molecules  within  the  interface.  In 
this  situation,  the  density  p is  the  variable  which  distinguishes  the  bulk  fluids  and  the 
intervening  interface.  In  this  role,  it  is  known  as  the  order  parameter.  In  contrast,  for  a 
binary  fluid  undergoing  spinodal  decomposition,  the  composition  c naturally  plays  the  role 
of  an  order  parameter  (Cahn  & Hilliard  1958).  Alternatively  it  is  possible  that  neither  the 
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density  or  the  composition  may  be  an  appropriate  or  convenient  order  parameter;  such  is 
the  case  in  solidification  models  (e.g.  Caginalp  1985,  1986,  Langer  1986).  Here  it  is  possible 
to  introduce  an  alternative  order  parameter,  the  so-called  phase  field  </>,  to  characterize  the 
phases.  The  phase-field  assumes  distinct  constant  values  in  each  bulk  phase  and  undergoes 
rapid  but  smooth  variation  in  the  interfacial  region.  The  phase  field  can  be  regarded  as 
a mathematical  device  that  allows  a reformulation  of  the  free-boundary  problem  and  has 
been  used  successfully  in  many  instances.  In  particular,  phase-field  models  of  solidification 
have  been  used  to  compute  complicated  realistic  interfacial  structures  such  as  those  present 
during  dendritic  growth  (Kobayashi  1993,  Wheeler  et  al  1993,  Warren  &;  Boettinger  1995) 
and  Ostwald  ripening  (Warren  &;  Murray  1996). 

Inherent  in  diffuse-interface  models  is  an  interfacial  width  which  is  characterized  by  the 
length  scale  over  which  the  order  parameter  changes.  By  considering  the  asymptotic  limit 
in  which  the  interfacial  width  is  small  compared  to  the  macroscopic  length  scale  associated 
with  the  motion  of  the  two  bulk  fluids  (i.e.  the  sharp-interface  limit)  the  diffuse-interface 
model  can  be  related  to  the  free-boundary  problem. 

In  Section  2 we  formulate  the  diffuse-interface  theory  of  a single-component  fluid  near  its 
critical  point.  Here  we  discuss  developments  for  equilibrium  and  nonequilibrium  situations 
and  also  review  the  applications  addressed  with  this  model.  In  Section  3 we  review  the 
developments  of  these  ideas  for  a binary  fluid.  In  Section  4 we  discuss  related  topics  including 
the  sharp-interface  limit  analyses,  computational  methods  for  fluid  interfaces  and  models  of 
miscible  fluids. 

2.  A SINGLE-COMPONENT  FLUID 

Diffuse-interface  models  of  a single- component  fluid  have  been  developed  largely  from  the 
perspective  of  critical  phenomena.  While  they  have  been  used  to  study  phenomena  associated 
with  the  critical  point,  they  have  also  been  applied  to  situations  away  from  the  critical 
point.  With  this  in  mind,  we  shall  review  in  this  section  the  ideas  and  applications  of 
diffuse-interface  models  of  a single-component  fluid. 
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Equilibrium 

We  begin  by  first  considering  the  equilibrium  state  of  a non-uniform  single- component  fluid. 
We  assume  that  an  isothermal  fluid  near  its  critical  point  has  associated  with  it  the  Helmholtz 
free  energy  functional  given  by 


Jy  [pKp,  T)  + ±K\VPf 


dV, 


(5) 


where  V is  a control  volume,  f(p,T)  is  the  bulk  free  energy  density  (per  unit  mass),  K is  a 
gradient  energy  coefficient  (assumed  for  simplicity  to  be  constant)  and  T is  the  temperature. 
In  a simple  model  the  term  pf(p,T)  is  assumed  to  take  the  form  of  a double-well  with 
respect  to  the  density  below  the  critical  temperature  and  a single-well  above  the  critical 
temperature.  The  square-gradient  term  is  associated  with  variations  of  the  density  and 
contributes  to  the  free  energy  excess  of  the  interfacial  region,  which  defines  the  surface 
energy  (Cahn  Sz  Hilliard  1958).  The  form  of  the  energy  density  can  be  interpreted  in  the 
context  of  statistical  mechanics  in  which  the  square  gradient  term  arises  from  attractive 
long-ranged  interactions  between  the  molecules  of  the  fluid  and  in  which  the  gradient  energy 
coefficient,  FT,  can  be  related  to  the  pair  correlation  function,  (e.g.  see  Irving  Sz  Kirkwood 
1950,  Bearman  Sz  Kirkwood  1958,  Yang  et  al  1976,  Bongiorno  et  al  1976,  Abraham  1979,  de 
Sobrino  Sz  Peternelj  1982  and  Davis  Sz  Scriven  1982  for  further  details). 

The  equilibrium  conditions  are  obtained  by  minimizing  T subject  to  a constraint  of 
constant  mass,  Ad,  where 


M = f pdV. 
Jv 

This  leads  to  the  Euler-Lagrange  equation 


(6) 


KV2p  - ( pf ),  + A = 0, 


(7) 


where  A is  the  Lagrange  multiplier  associated  with  conservation  of  mass.  We  observe  that 
the  integrand  of  T as  well  as  that  of  the  mass  constraint  are  independent  of  the  spatial 
coordinates.  Consequently,  it  follows  from  Noether’s  Theorem  (Goldstein  1980)  that  there 
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is  a corresponding  conservation  law  given  by 


V • T = 0, 


(8) 


where  T is  a second-rank  tensor  given  by 

T = CI-Vp 


dC 


5(Vp) 


(9) 


and  C = pf(p,T)  + \K\S7  p\2  — A p.  Using  the  Euler-Lagrange  equation  (7)  to  eliminate  the 
Lagrange  multiplier  we  find  that 


-p  + KpV2p  + -K\Vp\2 


I — KV p 0 Vp, 


(10) 


where  p = p2fp  has  been  identified  as  the  thermodynamic  pressure,  see,  e.g.,  Callen  (1985). 
Using  the  divergence  theorem,  equation  (8)  may  be  expressed  as 


J T ■ hdA  = 0,  (11) 

where  S is  the  boundary  of  a control  volume  with  unit  normal  vector  h.  This  is  the  gen- 
eralization to  three  dimensions  of  the  first  integral  of  the  Euler-Lagrange  equation  (7)  in 
one  dimension.  Equations  (8)  and  (11)  suggest  that  T represents  a stress  tensor  (up  to  an 
additive  divergence-free  contribution).  We  show  in  the  next  section  that  T represents  the 
reversible  part  of  the  stress  tensor.  Similar  equilibrium  conditions  have  been  obtained  by 
Blinowski  (1973a,  1973b)  from  the  point  of  view  of  elastic  fluids.  A review  article  covering 
a variety  of  aspects  of  this  and  related  theories  can  be  found  in  Davis  & Scriven  (1982).  It 
was  noted  by  Dunn  Sz  Serrin  (1985)  that  consistency  with  nonequilibrium  thermodynamics 
requires  a more  specific  form  for  the  capillary  tensor  than  that  used  originally  by  Korteweg 
(1901)  and  later  in  the  mechanical  equilibrium  theories  of  Aifantis  & Serrin  (1983a,  1983b). 

The  equilibrium  density  profile  p(z)  obtained  using  a van  der  Waals  equation  of  state 
(Callen  1985,  Rowlinson  and  Widom  1989)  represents  a smooth  transition  from  one  bulk 
density  to  the  other  (in  the  two  phase  region)  over  a length  scale  associated  with  the  gradient 
energy  coefficient.  Such  an  interface  has  associated  with  it  a surface  energy  given  by 
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Non- equilibrium 


We  now  pursue  the  non-equilibrium  situation  by  outlining  a thermodynamic  procedure  in- 
volving local  balances  of  mass,  linear  momentum,  energy  and  entropy  consistent  with  the 
inclusion  of  a square-gradient  energy  term  in  the  internal  energy  functional.  The  total  mass, 
M,  total  momentum,  V,  total  internal  energy,  S,  and  total  entropy,  S associated  with  a 


) are 

M = 

/ pdV, 

Jm 

(13a) 

V = 

f pvdV, 

Vn(t) 

(13b) 

e = 

(13c) 

5 = 

LrJV‘ 

(13d) 

where  v is  the  fluid  velocity,  e and  s are  the  internal  energy  and  entropy  per  unit  mass,  and 
Ke  is  the  gradient  (internal)  energy  coefficient  which  we  assume  to  be  constant.  Here  we 
have  for  simplicity  neglected  body  forces  such  as  gravity.  The  associated  physical  balance 
laws  can  be  expressed  as 


dM 

dt 

dV 


= 0, 


ay  _ r 
dt  J6 
dS 


5n(t) 


m ■ ndA , 


dt 

dS 


Ln(t)  ^ 


v • m ■ n — qE  • n]  dA, 


+ / qs-ndA=  [ 

dt  Jtn(t)  J n 


^oddV  > 0, 


(14a) 

(14b) 

(14c) 

(14d) 


rsn(t)  Jn(t) 

where  £f2(t)  is  boundary  of  0(f),  m is  the  stress  tensor,  qs  and  qs  are  the  internal  energy  and 
entropy  fluxes,  respectively,  and  is  the  volumetric  entropy  production.  The  quantities 
m,  qs  and  qs,  which  in  general  include  both  classical  and  nonclassical  contributions,  shall 
be  specified  below  such  that  their  forms  guarantee  that  i prod  is  non-negative,  as  required 
by  the  Second  Law  of  Thermodynamics.  Equation  (14a)  simply  represents  conservation  of 
mass.  Equation  (14b)  states  that  the  change  in  total  momentum  is  related  to  the  forces 
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on  the  boundary  (again  note  that  we  have  neglected  body  forces).  Equation  (14c)  states 
that  the  change  in  total  energy  is  related  to  the  rate  of  working  done  by  the  forces  on  the 
boundary  and  also  the  energy  flux  through  the  boundary.  Equation  (14d)  states  that  the 
change  in  total  entropy  plus  the  entropy  flux  through  the  boundary  must  be  equal  to  the 
entropy  production.  The  definitions  (13)  and  balance  laws  (14)  may  be  manipulated  to  show 
that 


\prod 


(m  — T)  : Vu 


+ 


(a  + &fv,)v(i) 


+ V - I qs  - - 


^V, 

T Dt  H 


(15) 


where  T is  the  reversible  part  of  the  stress  tensor  given  by  equation  (10)  with  K replaced  by 
Ke , and  we  have  used  the  thermodynamic  relationship  de  = Tds  + (p/ p2)dp.  The  following 
specifications  ensure  that  the  entropy  production  is  positive 


m 

= T + t, 

(16a) 

cC 

> 

cl|  to 

1 

E-i 

> 

1 

II 

(16b) 

qs 

i 

1 

II 

(16c) 

where  k is  the  thermal  conductivity  and  r is  the  viscous  stress  tensor  given  in  the  standard 
manner  as  r = 7/(V  *u)  + /z(Vu  + Vu7),  where  rj  and  p are  coefficients  of  viscosity  (e.g.  Batch- 
elor 1967).  This  prescription  for  m is  similar  to  that  postulated  by  Korteweg  (1901).  In 
this  formulation  it  is  assumed  that  the  viscosity  and  thermal  conductivity  are,  in  general, 
functions  of  the  density.  We  observe  that  the  energy  flux  qs  involves  both  the  classical 
contribution  corresponding  to  the  Fourier  Law  for  heat  conduction  (Carslaw  & Jaeger  1959, 
Kittel  Sz  Kroemer  1980)  and  a nonclassical  contribution.  This  nonclassical  contribution  was 
referred  to  as  ‘interstitial  working3  by  Dunn  & Serrin  (1985),  who  also  noted  that  there  is 
no  corresponding  nonclassical  term  in  the  entropy  flux.  When  a square-gradient  term  is  in- 
cluded in  the  definition  for  total  entropy  (13d),  a nonclassical  entropy  flux  arises  (e.g.  Wang 
et  al  1993,  Wheeler  et  al  1996,  Anderson  & McFadden  1996).  Using  the  above  forms  (16) 
for  the  stress  tensor  and  fluxes  the  local  balance  laws  may  then  be  written  as 


= -pV  • v, 


Dp 

Dt 
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= V ra,  (17b) 

= V-(fcVr)  + (-j>Z  + r)  : Vi;,  (17c) 

= V • (jfeVT)  + r : Vu,  (17d) 

These  equations  describe  viscous,  compressible,  nonisothermal  flow.  In  order  to  solve  these 
it  is  necessary  to  supply  an  equation  of  state.  In  the  isothermal  situation,  for  example,  one 
might  specify  the  pressure  p through  a van  der  Waals  equation  of  state. 

This  and  similar  models  have  been  developed  and  studied  by  a number  of  authors  from 
a variety  of  perspectives.  Fixman  (1967)  developed  a diffuse- interface  hydrodynamic  model 
in  which  the  reversible  part  of  the  stress  tensor  was  identified  using  a mechanical  principle. 
In  the  model  of  Felderhof  (1970),  a more  general  form  of  the  stress  tensor  was  obtained 
which  was  compatible  with  the  underlying  Lagrangian.  Langer  & Turski  (1973)  derived 
a diffuse-interface  model  that  they  showed,  through  a ‘coarse-graining’  argument,  could 
be  related  to  a molecular  model.  The  model  employed  by  Jasnow  &;  Vinals  (1996)  was 
derived  using  a Hamiltonian  description.  Jacqmin  (1996)  described  a model  which  also 
included  a wall  potential  to  model  the  interaction  between  the  fluid  and  a solid  boundary. 
Truskinovsky  (1993)  derived  a similar  model  which  also  included  an  additional  non-conserved 
order  parameter  and  its  gradients.  Antanovskii  (1996)  presented  a derivation  of  the  model 
based  on  a maximum  entropy  principle.  The  derivation  outlined  above  is  most  similar  to 
that  described  by  de  Sobrino  (1976),  Dunn  &:  Serrin  (1985),  Dunn  (1986)  and  Anderson 
Sz  McFadden  (1996,  1997).  The  work  of  de  Sobrino  begins  in  a more  general  framework 
and  invokes  symmetry  and  invariance  principles  to  simplify  the  gradient  dependence  of  the 
stress  tensor  to  that  which  is  described  above.  Dunn  & Serrin’s  approach  is  similar,  but 
mathematically  more  rigorous  and  is  given  from  the  viewpoint  of  rational  mechanics. 

A detailed  numerical  analysis  of  a simplified  version  of  the  diffuse-interface  model  has 
been  performed  by  Affouf  and  Caflisch  (1991),  who  present  numerical  solutions  representing 
phase  transitions,  shocks,  and  rarefaction  waves  connecting  far-field  states  and  analyze  their 
stability. 
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Applications 

CRITICAL  POINT  SCALING  LAWS  Extensive  analysis  using  renormalization-group 
techniques  have  been  performed  on  a diffuse-interface  model,  commonly  known  in  the  liter- 
ature as  ‘Model  H’  (Hohenberg  Sz  Halperin  1977),  which  describes  the  dynamics  of  a binary 
fluid  phase  transition  as  well  as  a single-component  fluid  near  its  critical  point  (e.g.  Halperin 
et  al  1974,  Siggia  et  al  1976,  and  Hohenberg  Sz  Halperin  1977).  Such  analyses  have  identified 
divergent  transport  coefficients  and  scaling  relations  associated  with  near-critical  fluids. 

SHEAR  FLOWS  IN  NEAR-CRITICAL  FLUIDS  Onuki  Sz  Kawasaki  (1979)  and  Onuki 
et  al  (1981)  have  studied  the  dynamics  of  a near-critical  fluid  in  a shear  flow  using  a model 
developed  by  Kawasaki  (1970).  They  investigated  the  regime  in  which  the  equilibrium  cor- 
relation length  (i.e.  the  interfacial  thickness)  exceeds  the  length  scale  associated  with  the 
shear  flow.  Among  their  findings  is  that  the  critical  fluctuations  of  classical  fluids  can  be 
drastically  altered  (e.g.  they  can  become  highly  anisotropic)  by  shear  flows.  These  ideas 
have  also  been  applied  to  polymers  under  shear  flows  (Helfand  Sz  Fredrickson  1989,  Onuki 
1989). 

CAPILLARY  WAVES  The  diffuse- interface  model  has  been  used  to  study  capillary 
waves  (Felderhof  1970,  Turski  Sz  Langer  1980,  de  Sobrino  Sz  Peternelj  1985  and  de  Sobrino 
1985).  These  authors  begin  with  a diffuse  planar  interface  in  equilibrium.  They  discuss 
capillary  waves  by  means  of  a linearized  theory  about  the  equilibrium  state.  The  linearized 
governing  equations  are  examined  in  the  long  wavelength  limit  and  the  dispersion  relation  for 
capillary  waves  obtained  from  the  classical  free-boundary  problem  (e.g.  Landau  Sz  Lifshitz 
1959)  is  recovered.  These  authors  recognized  the  importance  of  an  isentropic  treatment 
to  the  existence  of  capillary  waves  (in  the  context  of  a diffuse-interface  model)  although 
Felderhof  derived  the  result  in  the  isothermal  case  as  well  under  somewhat  more  restrictive 
conditions. 

MOVING  CONTACT  LINES  Seppecher  (1996)  established  the  governing  equations 
for  an  isothermal  viscous  flow  near  a moving  contact  line  on  a planar  solid  wall  using  a 
diffuse-interface  model.  In  particular,  a uniform  distribution  of  ‘double  forces5  was  used  to 
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describe  the  interaction  between  the  fluid  and  the  wall.  In  this  analysis,  the  contact-line 
problem  was  separated  into  three  regions:  an  external  (outer)  region  far  from  the  contact 
line  where  the  classical  theory  of  capillarity  applies;  an  inner  region  very  near  the  contact  line 
whose  dimensions  are  so  small  that  the  thickness  of  the  interface  cannot  be  neglected;  and  an 
intermediate  region  which  matches  the  inner  and  outer  regions.  The  flow  in  the  inner  region 
is  compressible  while  the  flows  in  the  intermediate  and  outer  regions  are  incompressible. 
Examples  of  the  density  field  and  flow  field  in  the  inner  region  are  shown  in  Figs.  1 and  2.  A 
key  result  of  this  paper  was  that  the  force  singularity  present  in  classical  continuum  models  of 
moving  contact  lines  (e.g.  Huh  & Scriven  1971,  Dussan  V.  & Davis  1974,  Dussan  V.  1979)  is 
not  present  when  the  interface  is  modeled  as  diffuse.  This  is  attributed  to  mass  transfer  across 
the  interfacial  region.  Additionally,  numerical  computations  reveal  a roughly  linear  increase 
in  the  apparent  contact  angle  with  the  contact-line  velocity  (no  contact-line  hysteresis  was 
considered),  which  agrees  with  the  general  trends  observed  experimentally  (see  Dussan  V. 
1979).  The  moving  contact-line  problem  has  also  recently  been  treated  computationally  by 
Jacqmin  (1996)  and  is  described  in  the  next  section. 

INTERNAL  WAVES  IN  A NEAR  CRITICAL  FLUID  Anderson  & McFadden  (1997) 
have  recently  employed  the  diffuse-interface  model  to  describe  internal  waves  in  a near- 
critical  fluid.  These  internal  waves  are  present  in  small  (centimeter-scale)  containers  owing 
to  the  large  compressibility  of  the  fluid  near  the  critical  point  and  have  been  observed 
experimentally  in  near-critical  Xenon  by  Berg  et  al  (1996).  In  the  experimental  work  of  Berg 
et  al,  internal  gravity  wave  frequencies  were  measured  both  above  the  critical  temperature 
where  a single  phase  exists  and  below  the  critical  temperature  where  two  phases,  separated 
by  an  interface,  exist.  The  theoretical  development  of  Berg  et  al  consisted  of  two  separate 
models;  one  above  the  critical  temperature  and  another  below  the  critical  temperature.  In 
contrast  to  these  classical  hydrodynamic  models,  the  diffuse-interface  approach  employed  by 
Anderson  &;  McFadden  allows  for  a single  model  to  be  applied  both  above  and  below  the 
critical  temperature.  In  their  diffuse-interface  model,  they  use  a van  der  Waals  equation 
of  state  to  obtain  a base-state  density  profile.  Upon  this  base-state  they  introduce  linear 
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perturbations  to  calculate  internal  wave  frequencies  for  temperatures  both  above  and  below 
the  critical  temperature.  Predictions  of  the  internal  wave  frequency  from  the  diffuse-interface 
model  are  compared  favorably  with  experimental  data  and  theoretical  results  of  Berg  et  al. 

DROPLETS  AND  NUCLEATION  Several  authors  have  investigated  the  nucleation 
of  droplets  (Blinowski  1974  and  DellTsola  et  al  1995,  1996).  The  focus  of  this  work  was  to 
ascertain  the  effects  of  interfacial  thickness  on  the  nucleation  conditions  of  a droplet  since,  at 
nucleation,  the  droplet  radius  may  be  comparable  to  the  interfacial  thickness.  Under  these 
circumstances,  the  classical  Laplace-Gibbs  theory  for  the  equilibrium  radius  of  a droplet  is 
called  into  question.  DellTsola  et  al  used  an  equilibrium  formulation  of  a diffuse-interface 
model  developed  in  earlier  work  (DellTsola  & Kosinski  1993)  to  study  nucleation  of  spherical 
droplets.  In  particular,  they  noted  that  for  microscopic  droplets  the  difference  in  mechanical 
pressures  inside  and  outside  the  droplet  (which  is  what  is  measured  experimentally)  is  not 
the  same  as  the  difference  between  the  thermodynamic  pressures  inside  and  outside  the 
droplet.  The  mechanical  pressure  involves  a (stress)  contribution  from  the  spatial  density 
variation  (e.g.  the  term  in  (10)  involving  p^72p)  which  at  the  center  of  the  microscopic  bubble 
is  important.  Based  on  these  ideas,  they  carefully  define  quantities  such  as  surface  tension 
and  droplet  radius  in  a way  which  generalize  the  classical  notions.  Relationships  between 
surface  tension,  droplet  radius  and  the  critical  nucleation  radius  are  obtained  using  a number 
of  different  equations  of  state.  Results  for  the  minimal  nucleation  radius  are  compared  with 
experimental  measurements. 

INSTABILITIES  OF  PLANAR  JETS  Nadiga  & Zaleski  (1996)  studied  in  the  instability 
of  a planar  jet  of  viscous,  compressible,  isothermal  liquid  issuing  into  its  surrounding  gas 
phase.  They  use  a van  der  Waals  equation  of  state  to  characterize  the  system.  Their 
calculations  focus  on  the  high  Reynolds  number  regime  (Re  = 800)  and  they  investigate  the 
effect  of  surface  tension  on  the  stabilization  of  the  jet. 

SPINODAL  DECOMPOSITION  IN  A PURE  FLUID  Nadiga  & Zaleski  (1996)  also 
considered  the  spinodal  decomposition  of  a single- component  fluid  rapidly  quenched  from  a 
temperature  above  its  critical  point  to  a temperature  below  the  critical  point.  These  compu- 
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tations  are  isothermal  and  two-dimensional.  They  find  that  the  inclusion  of  hydrodynamic 
effects  in  the  model  leads  to  a predicted  growth  rate  which  is  slightly  enhanced  from  the  case 
when  the  growth  is  limited  by  diffusion.  Numerical  computations  show  the  domain  growth 
and  plots  of  domain  size  versus  time  are  presented. 

3.  A BINARY  FLUID 


Non- equilibrium 

We  now  consider  the  situation  of  a binary  fluid  which  consists  of  two  components  A and 
B.  We  denote  the  composition  of  component  A,  expressed  as  a mass  fraction,  by  c.  In  this 
setting  the  composition  plays  the  role  of  an  order  parameter  which  distinguishes  the  different 
phases  of  the  fluid  and  in  this  way  is  analogous  to  the  density  in  the  single- component  fluid 
models.  As  we  shall  see  shortly,  however,  there  are  a number  of  subtle  differences  between 
the  single-component  and  binary  fluid  models. 

The  governing  equations  for  the  flow  of  a nonisothermal,  compressible  binary  fluid  may 
be  developed  in  a similar  manner  to  those  for  the  single- component  case  discussed  above  but 
involve  square-gradient  contributions  in  the  internal  energy  functional  from  the  composition 
rather  than  from  the  density.  This  procedure  gives  that  the  entropy  production  is 


(m-Tc):Vv  / _ Dc„  \ / 1 \ _ / 

5 = h ( <1E  + Ke~^S7c)  ‘ V \ ~ 4°  ' V ( 


+ V 


qE 


KeDc , 


9S---— — Vc  + ^-c 


where  Tc  is  the  reversible  part  of  the  stress  tensor  given  by 

To  = (-P  + ):Ke\ Vcl2)  I - Kb^c  ® Vc, 


(18) 


(19) 


qc  is  the  mass  flux  of  component  A and  /zc  is  the  generalized  chemical  potential  given  by 

(20) 


\ic  — t: V c. 


de 

dc 


The  specification  of  m,  qs,  qs  and  qc  to  ensure  positive  entropy  production  follows  in 
an  analogous  way  to  the  single- component  case.  In  particular,  qc  = -RV(/ic/T),  where 
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(21) 


D is  positive  and  represents  the  diffusion  coefficient.  Consequently,  the  equation  for  the 
composition  is  given  by 

9 Dt  \ T \dc  P 

This  is  well-known  Cahn-Hilliard  equation,  used  to  model  spinodal  decomposition  by  Cahn 
(1961)  and  Hilliard  (1970),  modified  to  account  for  fluid  motion.  The  equations  governing 
density,  velocity  and  temperature  (energy)  are  similar  to  those  for  a single- component  fluid 
given  above  by  equations  (17).  Two  subtle  differences  between  the  binary  fluid  equations 
and  the  single- component  fluid  equations  are  worth  mentioning.  First,  the  reversible  part  of 
the  stress  tensor  Tc  does  not  contain  the  counterpart  of  the  Laplacian  term  appearing  in  T, 
see  equation  (10).  This  is  because  the  order  parameter  c is  given  per  unit  mass  in  contrast 
to  the  density  p.  Second,  the  order  parameter  c is  governed  by  the  modified  Cahn-Hilliard 
equation  (21),  whereas  the  density  p is  governed  by  the  continuity  equation  (17a). 

Details  of  the  derivation  of  the  governing  equations  for  a binary  fluid  have  been  dis- 
cussed by  a number  of  authors.  Blinowski  (1975)  considered  the  binary  fluid  case  and  also 
more  general  multi- component  systems.  Starovoitov  (1994)  derived  a binary  fluid  diffuse- 
interface  model  using  a virtual  power  method.  Antanovskii  (1995)  derived  the  equations  for 
nonisothermal,  viscous,  quasi-compressible  flow.  His  derivation  was  based  on  a maximum 
entropy  principle,  similar  to  the  approach  described  here,  but  used  virtual  work  arguments 
to  identify  the  forms  of  the  stress  tensor  and  fluxes.  Gurtin  et  al  (1996)  derived  a model  for 
an  isothermal,  incompressible  flow  using  microforce  balance  laws.  Jasnow  & Vinals  (1996) 
illustrated  the  derivation  of  the  equations  for  a binary  fluid  using  a Hamiltonian  formalism. 
Although  their  derivation  was  given  for  an  isothermal,  inviscid,  incompressible  fluid,  the 
model  was  extended  to  the  case  where  the  fluid  was  viscous  and  the  temperature  field  varied 
slowly  in  time.  A similar  account  can  be  found  in  Jacqmin  (1996),  who  discussed  the  model 
in  a ‘potential  form’  as  well  as  making  note  of  its  relation  to  the  above  ‘stress  form’.  The 
merits  of  these  two  equivalent  formulations  were  discussed  in  terms  of  specific  applications. 
Lowengrub  &;  Truskinovsky  (1997)  presented  a thorough  derivation  of  the  diffuse-interface 
model  based  on  entropy  production  and  paid  particular  attention  to  the  difference  between 
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compressible  and  quasi-incompressible  fluids.  The  quasi-incompressible  situation  describes 
the  case  where  the  fluid  density  is  independent  of  the  pressure.  However,  the  bulk  states 
may  have  different  densities  and  the  flow  in  the  interfacial  region  is  in  general  non-solenoidal 
(V  -v  ^ 0),  resulting  in  an  expansion  or  contraction  flow  upon  phase  transformation.  These 
authors  argued  that  within  the  context  of  quasi-incompressibility,  the  appropriate  thermody- 
namic description  is  in  terms  of  a Gibbs  free  energy,  in  which  the  pressure  is  an  independent 
variable  determined  by  the  transport  equations  rather  than  a quantity  determined  thermo- 
dynamically. 

Applications 

THERMOCAPILLARY  FLOWS  Antanovskii  (1995)  used  his  nonisothermal  binary  fluid 
model  to  compute  one-dimensional  thermocapillary  flow  in  a gap.  In  this  situation,  two  fluid 
phases,  characterized  by  their  different  compositions,  were  separated  by  a planar  diffuse 
interface  along  which  a temperature  gradient  was  imposed.  Calculations  for  different  values 
of  interfacial  layer  thickness  and  viscosity  ratios  were  presented.  Jasnow  &:  Vinals  (1996) 
investigated  thermocapillary  migration  of  droplets  of  one  phase  in  the  surrounding  phase. 
They  focused  on  drops  with  radii  on  the  order  of  ten  times  the  correlation  length.  Their 
calculations  show  the  motion  of  the  droplet  through  the  temperature  gradient  as  a function 
of  time  and  the  dependence  of  its  velocity  on  this  temperature  gradient.  Also  shown  is  a 
sequence  in  which  two  droplets  coalesce. 

SPINODAL  DECOMPOSITION  Gurtin  et  al  (1996)  considered  spinodal  decomposition 
occurring  in  an  isothermal  binary  fluid.  Their  computations,  which  begin  with  an  initial 
random  distribution  of  the  composition,  show  the  coarsening  process  explicitly.  They  note 
that  the  main  effect  of  the  hydrodynamic  interactions  on  this  process  is  the  flow-induced 
coalescence  of  droplets.  Their  calculations  specifically  show  the  flow  associated  with  these 
coalescence  phenomena.  They  compute  the  domain  (structure)  size,  compare  with  classical 
predictions,  and  find  that  for  long  times  the  growth  is  faster  than  the  classical  scaling  law 
for  coarsening  by  purely  diffusive  mechanisms.  Jasnow  & Vinals  (1996)  consider  a similar 
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situation  but  include  the  effects  of  a nonuniform  temperature  field.  Their  results,  displayed 
in  the  time-sequence  in  Fig.  3,  show  the  spinodal  decomposition  of  a binary  fluid  in  a 
rectangular  cell  across  which  a vertical  temperature  gradient  is  imposed  (the  cell  is  hottest 
at  the  bottom).  Here,  an  initially  random  composition  distribution  (top  frame)  undergoes 
coarsening  and  domain  growth.  At  an  intermediate  time  (middle  frame)  and  a later  time 
(bottom  frame),  their  calculations  show  nonuniform  coarsening  wherein  the  smaller  (larger) 
scales  occur  in  the  warmer  (cooler)  regions. 

MIXING  AND  INTERFACIAL  STRETCHING  Chella  k Vinals  (1996)  compute  the 
mixing  and  interfacial  stretching  of  an  isothermal,  incompressible  binary  fluid  in  a shear 
flow  with  equal  densities  and  viscosities  in  the  two  phases.  An  initially  planar  interface  is 
shown  to  distort  under  an  imposed  shear  flow.  The  wrapping  up  of  the  two  phases  is  shown 
for  different  values  of  the  capillary  number  (surface  energy).  The  amount  of  interfacial 
stretching  is  computed  and  compared  with  an  analytical  solution  for  that  associated  with  a 
passive  scalar.  An  increase  in  the  surface  energy  corresponds  to  a decrease  in  the  amount 
of  stretching  and  in  this  way  opposes  the  effect  of  the  shear  flow.  Further  calculations  by 
Chella  and  Vinals  of  flow  in  a driven  cavity  are  shown  in  Fig.  4.  A number  of  breakup  and 
coalescence  events  can  be  seen  in  this  sequence. 

DROPLET  BREAKUP  Jacqmin  (1996)  calculated  the  breakup  of  an  inviscid  fluid  in 
the  context  of  a two-dimensional  isothermal  model.  The  initially  elongated  droplet  relaxes, 
oscillates  and  breaks  apart  into  two  separate  droplets.  In  these  calculations,  the  diffusion 
coefficient  ( D in  our  notation)  was  made  velocity  dependent.  In  addition  to  plots  showing 
the  actual  droplet  breakup  process,  the  kinetic  and  surface  energy  evolutions  are  also  given. 
In  the  case  where  the  flow  is  strongly  damped  by  using  a larger  value  of  D,  the  energy  of 
the  droplet  decreases  rapidly  and  breakup  does  not  occur. 

WAVE-BREAKING  AND  SLOSHING  Jacqmin  (1996)  also  applied  the  diffuse-interface 
model  to  a large- deformation  sloshing  flow  in  a two-dimensional  rectangular  domain.  Inter- 
action with  the  solid  wall  of  the  container  is  modeled  using  a wall  potential  associated  with 
a 90°  contact  angle  between  the  wall  and  the  interface.  The  two-phase  fluid  in  a uniform 
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(vertical)  gravitational  field  was  subjected  to  an  oscillating  horizontal  acceleration  whose 
amplitude  increased  linearly  in  time.  The  motion  of  the  initially  planar  interface  between 
the  two  stably  stratified  phases  is  computed  once  the  horizontal  accelerations  begin.  It  is 
shown  that  the  shape  of  interface  becomes  highly  nonlinear  and  several  breakup  and  coa- 
lescence events  can  be  seen  to  occur.  It  is  noted  that  the  prediction  of  coalescence  events 
may  occur  more  quickly  than  is  physically  realistic  since  a high  level  of  accuracy  is  needed 
to  resolve  the  draining  layer  between  the  two  coalescing  phases. 

MOVING  CONTACT  LINES  Jacqmin  (1996)  investigated  the  fluid  motion  near  a 
moving  contact  line.  These  steady  calculations  show  the  existence  of  a dynamic  contact 
angle  associated  with  the  interface  away  from  the  immediate  vicinity  of  the  contact  line. 
Also  observed  is  a streamline  (similar  to  that  observed  in  the  experiments  by  Dussan  V.  &; 
Davis  (1974))  issuing  from  the  contact  line  region  into  the  displaced  fluid. 

NUCLEATION  Lowengrub  & Truskinovsky  (1997)  considered  nucleation  (and  annihi- 
lation) of  an  isothermal,  spherically  symmetric  equilibrium  droplet.  In  particular,  they  used 
the  diffuse-interface  model  to  describe  the  situation  where  the  droplet  size  was  comparable 
to  the  interface  thickness.  In  the  incompressible  case,  where  the  density  was  uniform  ev- 
erywhere, they  performed  an  analysis  of  the  spherically  symmetric  Cahn-Hilliard  equation 
with  a free-energy  density  composed  of  piecewise  parabolas  to  find  new  analytic  solutions  for 
the  compositions.  In  the  compressible  case  they  proceeded  numerically  and  concluded  from 
their  results  that  compressibility  has  little  effect  on  the  interfacial  structure  of  the  droplet. 

4.  RELATED  TOPICS 

Sharp- Interface  Limit 

As  noted  in  Section  2,  the  diffuse-interface  models  may  be  applied  away  from  the  critical 
point,  where  the  interfacial  thickness  approaches  that  of  a ‘sharp’  boundary.  The  use  of 
the  diffuse-interface  models  in  these  regimes  may  be  justified  by  demonstrating  that  they 
approach  asymptotically  the  free-boundary  formulation.  This  is  achieved  by  adopting  what 
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is  commonly  known  as  the  ‘sharp-interface  limit’. 

There  has  been  a significant  effort  in  the  context  of  phase-field  models  of  solidification 
to  address  the  sharp-interface  limit  and  to  compare  the  model  to  well-known  results  of  the 
free-boundary  problem  (e.g.  Caginalp  1989,  Braun  et  al  1994). 

One  of  the  key  features  in  the  diffuse-interface  models  described  here,  which  is  not  present 
in  the  phase-field  solidification  models  studied  to  date,  is  the  (vector)  momentum  equation 
which  involves  the  distributed  capillary  stresses.  Consequently,  much  of  the  emphasis  in 
terms  of  the  sharp-interface  analyses  of  the  diffuse-interface  model  has  been  to  recover  the 
classical  interfacial  boundary  conditions  associated  with  the  stress  balance  at  the  interface. 
Here  we  describe  the  efforts  in  this  area  and  outline  a simple  reduction  of  the  momentum 
balance  to  the  interfacial  stress  jump  using  a pillbox  argument. 

Antanovskii  (1995)  derived  from  his  diffuse-interface  model  of  a binary  fluid  the  special 
cases  of  the  classical  hydrostatic  balance  for  a flat  interface  in  equilibrium  and  the  Young- 
Laplace  equation  for  a spherical  interface  in  equilibrium.  The  latter  case  has  also  been 
considered  for  a single- component  model  by  Blinowski  (1979).  Nadiga  and  Zaleski  (1996) 
also  confirmed  numerically  that  their  diffuse-interface  model  accurately  recovered  the  clas- 
sical results  for  a flat  interface  and  for  a liquid  droplet  in  equilibrium.  Jasnow  and  Vinals 
derived  from  the  capillary  term  in  their  momentum  equation  the  appropriate  sharp-interface 
tangential  and  normal  forces  when  the  surface  tension  was  a slowly  varying  function  along  the 
interface.  More  detailed  analytical  approaches  have  been  addressed  by  Starovoitov  (1994), 
Anderson  & McFadden  (1996)  and  Lowengrub  & Truskinovsky  (1997). 

To  illustrate  these  ideas,  we  apply  a pillbox  argument  to  the  momentum  balance  (17b)  to 
show  how  the  classical  stress  balance  at  a fluid-fluid  interface  can  be  derived  from  the  diffuse- 
interface  model.  We  define  a small  parameter  e measuring  the  thickness  of  the  interface  by 
Ke  = e2K.  We  then  consider  the  surface  5/  defined  by  the  contour  of  density  upon  which  the 
interfacial  region  collapses  in  the  limit  e — ► 0 and  define  a pillbox  enclosing  a portion  of  this 
surface  at  a fixed  point  in  time  in  such  a way  that  the  top  of  the  pillbox  is  above  the  surface 
at  a height  r = 6 and  the  bottom  of  the  pillbox  is  below  the  surface  at  a height  r = —8. 
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Here,  r is  a local  coordinate  normal  to  the  interface.  The  key  limit  in  the  pillbox  argument 
is  that  e < ^ < I where  L is  an  0(1)  length  scale  associated  with  the  outer  flow.  In  this 
limit,  the  volume  of  the  pillbox  becomes  negligible  on  the  outer  scales  but  the  variations  in 
the  density,  which  define  the  interfacial  region,  occur  over  a region  fully  contained  within  the 
pillbox.  Also  in  this  limit,  the  unit  normal  vectors  on  the  top  and  bottom  of  the  pillbox  are 
n and  — h respectively  while  the  unit  normal  on  the  side  is  given  by  m (note  that  rh  • h = 0). 
We  integrate  equation  (17b)  to  obtain 

0 = Jv  + v ■ (pv®  v)  - V • dV,  (22) 

where  we  have  used  (17a)  and  the  fact  that 

Dv  d(pv)  . _ _» 

P15i  = ~dt  + v ' ® v)  ■ (23) 


Next,  we  note  that 

/ ~dV  ~ f (pv)vi  ■ risdS , (24) 

J Vp  UL  J Sp 

where  vj  is  the  velocity  of  the  surface  Si  described  above  and  Sp  denotes  the  surface  of 
the  pillbox.  This  result  follows  by  translating  to  a frame  moving  with  the  interface  so  that 
d(pv)/dt  = d(pv)/dt'  — vi  ■ V(pu)  = d(pv)/dt'  — V • (pvi  <g>  v)  + pvV  ■ fi/.  The  terms 
d(pv)/dt'  + pvV  • vi  are  bounded  and  hence  do  not  contribute  as  the  pillbox  volume  goes  to 
zero.  We  use  this  and  the  divergence  theorem  to  obtain  from  equation  (22) 


0 = J \^pv(y  — vi)  • hs  — ns  • m]  dS.  (25) 

We  further  argue  that  the  fluid  velocity  terms  are  bounded  (so  that  they  do  not  contribute 
to  the  integral  over  the  side  surface  of  the  pillbox)  and  that  the  nonclassical  terms  do  not 
contribute  to  the  upper  and  lower  surfaces  of  the  pillbox  (since  e <C  6)  so  that 

0 = [ \pv(y  — vi)  • n — n • (—  pi  + t)|+cL4  — f rh-TdS , (26) 

where  A is  the  portion  of  5/  within  the  pillbox.  Now,  local  to  the  interfacial  region  we  have 
Vp  ~ pTh  and  V2p  ~ prr  to  leading  order.  Then  m • T ~ (m  • T • m)rh  so  that 

0 = J \j)v(y  — vi)  ■ n — n • {—pi  + t)]  dA  — j)  J (m  • T ■ rh)rhdrdl , (27) 


-19- 


DIFFUSE-INTERFACE  METHODS  IN  FLUID  MECHANICS 


April  1995 


where  C is  the  contour  defined  by  the  intersection  of  5/  and  the  pillbox  surface  and  dl  is  the 
increment  of  arclength  along  C . We  next  define  the  scalar 


7 = 


/OO 

(m  • T • rh)di 

-OO 


(28) 


which  can  be  shown  to  be  equal  to  the  excess  Kramer’s  potential,  that  is,  the  surface  energy 
(Anderson  & McFadden  1996).  We  then  apply  the  surface  divergence  theorem  (Weatherburn 
1925)  which  allows  us  to  write 


j)  'yrridl  = J [V57  — 7 (Vs  • h)  n]  dA, 


(29) 


where  V57  is  the  surface  gradient  of  7 and  Vs -7 i,  the  surface  divergence  of  n,  is  the  interfacial 
curvature,  JC.  This  gives 


0 = J | \pv(y  — vj)  ■ n — n • (—pi  + r)J + — [V57  — 7 ATn]  j dA. 


(30) 


Finally,  noting  that  the  area  of  integration  is  arbitrary  yields 


(-pI+T)]^ 


V57  - 7 /Cn, 


(31) 


which  is  the  classical  stress  balance  at  a fluid-fluid  interface  (Delhaye  1974).  Note  that  the 
first  term  on  the  left-hand  side  represents  a jump  in  the  momentum  across  the  interface  and 
is  zero  when  there  is  no  mass  flux  across  the  interface. 


Other  Computational  Methods 

The  diffuse- interface  models  share  common  features  with  a number  of  methods  developed 
from  a more  computational  point  of  view.  We  outline  several  of  these  methods  below  and 
highlight  some  of  their  applications. 

The  volume  of  fluid  (VOF)  method  (Hirt  & Nichols  1981,  Hyman  1984,  Tsai  & Yue 
1996)  is  a numerical  approach  to  the  free-boundary  problem  in  which  an  auxiliary  function 
F,  which  distinguishes  one  fluid  from  another,  is  introduced  in  order  to  identify  the  shape 
and  evolution  of  the  free  boundary.  This  function  satisfies  DF / Dt  = 0;  that  is,  it  is  advected 
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with  the  flow.  Each  computational  cell  has  associated  with  it  a value  of  F and  those  cells 
which  take  on  a value  between  two  bulk  values  of  F are  assumed  to  contain  part  of  the 
interface.  The  normal  to  the  interface  in  the  cell  is  determined  by  the  direction  of  the 
largest  local  gradient  and  the  position  of  the  interface  in  the  cell  is  arranged  so  that  F is  the 
fractional  volume  of  fluid  in  the  cell.  The  VOF  method,  like  the  diffuse- interface  approach, 
allows  a straightforward  description  of  flows  involving  complicated  boundary  shapes  and 
topological  changes.  In  contrast  to  the  diffuse-interface  approach,  however,  free  surface 
boundary  conditions  in  the  VOF  method  must  still  be  applied  at  the  free  boundary. 

Brackbill  et  al  (1992)  developed  a ‘continuum  surface  force5  model  wherein  they  iden- 
tified a volume  force  which  represents  surface  tension  spread  over  a small  but  finite  three- 
dimensional  interfacial  domain.  This  volume  force  was  related  to  a ‘color5  function  which,  for 
example,  can  represent  density  for  incompressible  flows.  The  defining  characteristics  of  this 
volume  force  are  that  it  gives  the  correct  surface  force  in  the  limit  of  a sharp  interface  and 
is  nonzero  only  in  the  interfacial  region.  This  approach  allows  a single-domain  description 
of  the  two-fluid  system  and  does  not  require  the  direct  application  of  boundary  conditions, 
which  are  built  into  the  governing  equations.  A number  of  numerical  results  are  presented 
for  both  static  and  dynamic  situations:  (z)  a static  drop  (rod),  (ii)  a nonequilibrium  rod 
upon  which  capillary  waves  move  along  the  surface,  (m)  the  Rayleigh-Taylor  instability  and 
the  associated  interfacial  deformation,  ( iv ) flow  induced  by  wall-adhesion  whereby  the  fluid 
conforms  to  an  imposed  equilibrium  contact  angle  on  the  wall,  and  (n)  jet-induced  tank 
mixing  and  liquid  reorientation  in  microgravity  environments. 

Another  related  approach  has  been  developed  by  Unverdi  and  Tryggvason  (1992a,  1992b). 
Their  approach  is  a front-tracking  technique  which  employs  a numerically- diffuse  description 
of  the  interface.  They  construct  an  indicator  function,  based  on  the  known  position  of  the 
(sharp)  interface,  which  identifies  fluid  properties  such  as  density  and  viscosity.  This  function 
is  then  artificially  spread  out  over  a small  region  on  the  scale  of  the  computational  mesh  size, 
allowing  the  fluid  properties  to  vary  smoothly  through  this  interfacial  region.  The  surface 
force  (i.e.  surface  tension)  is  also  distributed  over  this  interfacial  region  so  that  a single- 
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domain  approach  can  be  used  to  calculate  the  flow.  This  flow  then  determines  how  the 
interface  is  advected.  In  Unverdi  and  Tryggvason  (1992a),  both  two-  and  three-dimensional 
multiple  bubble  motion  and  interaction  are  presented.  In  Unverdi  and  Tryggvason  (1992b), 
the  Rayleigh-Taylor  instability  in  two  and  three  dimensions  and  also  bubble-bubble  interac- 
tion in  three  dimensions  are  computed.  Nobari  et  al  (1996)  have  recently  used  this  approach 
to  compute  head-on  collision  of  two  viscous  liquid  droplets  with  surface  tension.  Here,  rup- 
ture is  modeled  by  artificially  removing  the  thin  film  between  the  two  drops  at  a prescribed 
time.  They  found  that  if  no  rupture  takes  place,  the  drops  always  rebound,  but  that  when 
rupture  occurs  the  drops  may  later  split.  This  method  has  been  extended  by  Juric  and  Tryg- 
gvason (1995,  1996)  to  describe  flows  in  the  presence  of  phase  change.  Here,  they  applied 
the  model  to  vapor  bubble  dynamics  and  film  boiling. 

Another  highly  successful  computational  scheme  that  has  been  applied  to  interfacial 
motion  is  the  level  set  method,  see  Osher  & Sethian  (1988)  and  Sethian  (1996).  With  this 
method,  the  interface  is  represented  as  a level  set  of  a smooth  auxiliary  function  which  is 
computationally  analogous  to  the  order  parameter  used  in  diffuse-interface  descriptions.  An 
advantage  of  the  level  set  method  is  that  the  interface  remains  sharp  in  this  formulation, 
which  eliminates  the  need  for  added  numerical  resolution  in  the  direction  normal  to  the 
interface.  Within  the  context  of  fluid  mechanics  and  two-fluid  flows,  surface  tension,  for 
example,  is  represented  in  the  momentum  equations  as  a distributed  force  through  the  use 
of  a smoothed  delta-function  (Sussman  et  al  1994,  Chang  et  al  1996).  The  momentum 
equations  are  then  used  to  compute  the  flow  over  the  whole  domain  and  the  level  set  function 
is  advected  with  the  flow.  After  a normalization  procedure,  the  level  set  function  determines 
the  new  position  of  the  interface.  Mulder  et  al  (1992)  applied  the  level  set  approach  (without 
surface  tension)  to  study  the  Rayeigh-Taylor  and  Kelvin-Helmholtz  instabilities  within  the 
context  of  compressible  gas  dynamics.  Sussman  et  al  (1994)  used  this  approach  to  model 
incompressible  two-dimensional  rising  bubbles  and  falling  drops,  which  show  large  distortions 
and  topological  changes.  They  also  show  two  impacting  droplets  as  well  as  a single  droplet 
impacting  a surface.  Chang  et  al  (1996)  also  applied  a level  set  approach  for  incompressible 
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fluids  to  several  topologically  complex  flows.  They  presented  computations  of  two  merging 
fluid  bubbles  of  equal  density  and  also  two  merging  fluid  bubbles  of  different  density.  Further, 
they  investigated  the  Rayleigh-Taylor  instability  of  an  initially  motionless,  vertical  column 
of  fluid  and  show  calculations  of  the  subsequent  vortex  sheet  roll-up  phenomena. 

Miscible  Fluids 

The  ideas  that  are  involved  in  the  diffuse-interface  models  described  in  the  preceding  sec- 
tions are  similar  in  many  ways  to  those  for  miscible  fluids.  In  particular,  fluid  stresses 
(i.e.  Korteweg  stresses)  that  arise  due  to  concentration  and  density  gradients  at  the  interface 
between  two  miscible  fluids  lead  to  the  notion  of  surface  tension  between  miscible  fluids. 
Joseph  (1990)  has  investigated  these  ideas  both  experimentally  and  theoretically.  He  per- 
formed a number  of  experiments  and  highlighted  other  experimental  work  in  which  miscible 
liquid  droplets  rising  or  falling  in  another  liquid  exhibit  capillary-type  effects;  that  is,  their 
shapes  are  consistent  with  the  presence  of  surface  tension  on  the  interface.  In  the  theoretical 
development,  the  equations  governing  the  motion  of  the  fluid  are  similar  to  those  presented 
in  the  above  section  on  binary  fluids.  They  are  the  continuity  equation,  the  Navier-Stokes 
equations  modified  to  account  for  the  gradient  stresses,  the  heat  (energy)  equation,  and  a 
standard  diffusion  equation,  rather  than  the  modified  Cahn-Hilliard  equation  (21),  that  de- 
scribes the  evolution  of  the  composition.  Besides  the  difference  in  the  equations  governing 
the  composition,  another  key  difference  between  the  description  given  by  Joseph  (1990)  and 
that  for  a binary  fluid  described  in  section  3 is  that  his  model  described  a generalized  (or 
quasi)  incompressible  fluid.  Here  it  is  assumed  that  the  density  is  a function  of  compo- 
sition and  temperature  but  is  independent  of  the  pressure.  Although  the  density  then  is 
unchanged  by  pressure  variations,  the  associated  flow  may  still  be  nonsolenoidal  (V  • v ^ 0) 
due  to  composition  or  temperature  variations.  This  notion  has  been  carefully  adopted  into 
the  diffuse-interface  model  for  a binary  fluid  by  Lowengrub  & Truskinovsky  (1997). 

The  theoretical  model  described  by  Joseph  (1990)  has  been  analyzed  in  more  detail  by 
Galdi  et  al  (1991).  These  authors  rework  the  equations  and  identify  a new  solenoidal  velocity 
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which  is  a linear  combination  of  the  original  velocity  field  and  the  concentration  gradient. 
Associated  with  this  new  velocity  is  a pressure  field  that  is  a linear  combination  of  the 
original  pressure  and  the  divergence  of  the  original  velocity.  They  address  the  linear  and 
energy  stability  of  a quiescent,  vertically  (unstably)  stratified  incompressible  fluid  in  which 
Korteweg  stresses  arise  due  to  composition  gradients.  This  is  the  analog  of  the  classical 
Benard  problem  with  the  exception  that  the  authors  do  not  immediately  invoke  the  Boussi- 
nesq  approximation.  They  find  that  the  stability  results  depend  strongly  on  the  value  of  the 
coefficient  in  the  Korteweg  stress  (and  in  particular  on  its  sign).  When  its  sign  is  chosen 
consistent  with  that  in  the  capillary  stress  term  described  in  Section  2,  and  when  its  effect 
is  strong  enough,  an  unconditionally  stable  base  state  is  predicted. 

An  notable  idea  described  in  the  paper  by  Joseph  (1990)  is  that  of  a dynamic  (or  nonequi- 
librium) surface  tension  between  two  miscible  fluids.  That  is,  although  for  two  miscible  fluids 
it  is  not  clear  that  one  can  define  an  equilibrium  surface  tension,  based  on  the  idea  of  stresses 
associated  with  gradients  in  density  or  composition,  transient  or  dynamic  surface  tension  be- 
tween two  mixing  phases  can  be  studied.  Joseph  et  al  (1996)  consider  this  situation  in  detail. 
In  this  paper,  they  carry  out  an  analysis  using  the  model  of  Joseph  (1990)  of  transient  or 
dynamic  interfacial  tension  during  the  smoothing  of  an  initial  discontinuity  of  composition 
across  plane  and  spherical  surfaces  separating  two  miscible  liquids.  It  is  found  that  the 
dynamic  interfacial  tension  decays  in  time  like  t-1/2  and  has  contributions  from  Korteweg 
stresses  and  from  the  expansion  velocity  (i.e.  the  nonsolenoidal  part)  which  also  involves 
the  rate  of  change  of  viscosity  with  composition.  For  a plane  mixing  front,  diffusion  has 
a similarity  solution  and  they  show  that  there  is  no  associated  pressure  jump  due  to  the 
Korteweg  terms  (i.e.  surface  tension  does  not  lead  to  a jump  in  pressure  when  the  curvature 
of  the  interface  is  zero).  The  only  pressure  difference  across  the  mixing  front  is  due  to  the 
hydrostatic  pressure  difference  (i.e.  due  to  gravity)  across  the  layer  since  its  thickness  is 
growing  with  time.  For  an  initially  spherical  droplet  placed  in  another  (miscible)  fluid  they 
identify  a pressure  jump  associated  with  two  effects:  the  first  is  due  to  the  Korteweg  stress 
and  the  second  is  due  to  the  expansion  velocity  and  is  proportional  to  the  rate  of  change  of 
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viscosity  with  composition.  The  latter  term  can  take  on  either  sign,  depending  on  how  the 
viscosity  depends  on  the  composition. 

5.  SUMMARY 

We  have  reviewed  the  development  and  application  of  diffuse-interface  models  for  both  single- 
component and  binary  fluids.  These  models  have  foundations  in  statistical  mechanics,  ki- 
netic theory,  mechanical  theories,  and  nonequilibrium  thermodynamics.  They  provide  an 
alternative  approach  to  the  classical  hydrodynamic  free-boundary  problem  and  have  been 
used  successfully  in  many  applications.  They  should  continue  to  play  an  important  role,  in 
concert  with  other  theoretical  and  experimental  efforts,  in  the  understanding  of  hydrody- 
namic phenomena  associated  with  critical  fluids  and  other  flows  involving  complex  interface 
morphologies  and  topological  changes. 
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Figure  1:  Moving  contact  line:  This  figure  shows  the  constant  density  contours 
in  the  inner  region.  The  frame  of  reference  is  fixed  on  the  contact  line  region 
so  that  the  bottom  plate  moves  from  left  to  right  at  an  imposed  velocity.  In 
this  case  the  imposed  static  contact  angle  at  the  wall  is  approximately  54° 
while  the  dynamic  contact  angle  (away  from  the  immediate  vicinity  of  the 
contact  line)  is  approximately  125°.  The  parameter  values  consistent  with  the 
notation  given  in  Seppecher  (1996)  are  R = 20,  um  = 2,  K = 10,  g = —0.3, 
and  Ca  = 40  x 10~3.  (This  figure  was  provided  by  P.  Seppecher). 
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Figure  2:  Moving  contact  line:  This  figure  shows  the  streamlines  in  the  inner 
region  associated  with  the  density  contours  shown  in  Fig.  1.  The  frame  of 
reference  is  fixed  on  the  contact  line  region  so  that  the  bottom  plate  moves  from 
left  to  right  at  an  imposed  velocity.  (This  figure  was  provided  by  P.  Seppecher). 
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Figure  3:  Spinodal  decomposition:  This  time  sequence  (top  to  bottom)  shows 
spinodal  decomposition  of  a binary  fluid  in  a rectangular  cell  across  which  is 
imposed  a vertical  temperature  gradient  (the  cell  is  hottest  at  the  bottom). 
The  composition  is  indicated  in  greyscale.  These  computations  show  the  sys- 
tem evolving  from  an  initially  random  composition  distribution  in  the  top 
frame,  to  one  in  the  bottom  frame  where  the  domain  structure  size  varies  with 
the  temperature.  (This  figure  was  provided  by  D.  Jasnow  and  J.  Vinals). 


-38- 


DIFFUSE-INTERFACE  METHODS  IN  FLUID  MECHANICS 


April  1995 


Figure  4:  Flow  in  a driven  cavity:  This  time  sequence  (across  the  top,  left 
to  right  and  across  the  bottom,  left  to  right)  shows  the  order  parameter  (in 
greyscale)  representing  the  binary  fluid  composition.  The  flow  is  set  up  by 
an  imposed  velocity  (left  to  right)  on  the  bottom  surface  of  the  cavity.  The 
densities  and  viscosities  of  both  phases  are  equal.  (This  figure  was  provided 
by  R.  Chella  and  J.  Vinals). 
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